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Abstract. Results on the Prandtl-Blasius type kinetic and thermal boundary layer 
thicknesses in turbulent Rayleigh-Benard convection in a broad range of Prandtl 
numbers are presented. By solving the laminar Prandtl-Blasius boundary layer 
equations, we calculate the ratio of the thermal and kinetic boundary layer thicknesses, 

1/2 

which depends on the Prandtl number Vr only. It is approximated as 0.588'Pr ' 
for Vr < Vr* and as 0M2Vr~^'^ for Vr* < Vr, with Vr* = 0.046. Comparison 
of the Prandtl-Blasius velocity boundary layer thickness with that evaluated in the 
direct numerical simulations by Stevens, Verzicco, and Lohse (J. Fluid Mech. 643, 495 
(2010)) gives very good agreement. Based on the Prandtl-Blasius type considerations, 
we derive a lower-bound estimate for the minimum number of the computational 
mesh nodes, required to conduct accurate numerical simulations of moderately high 
(boundary layer dominated) turbulent Rayleigh-Benard convection, in the thermal and 
kinetic boundary layers close to bottom and top plates. It is shown that the number 
of required nodes within each boundary layer depends on Mi and Vr and grows with 
the Rayleigh number TZa not slower than ^ TZaP'^^. This estimate agrees excellently 
with empirical results, which were based on the convergence of the Nusselt number in 
numerical simulations. 
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1. Introduction 

Rayleigh-Benard (RB) convection is the classical system to study properties of thermal 
convection. In this system a layer of fluid confined between two horizontal plates 
is heated from below and cooled from above. Thermally driven flows are of utmost 
importance in industrial applications and in natural phenomena. Examples include the 
thermal convection in the atmosphere, the ocean, in buildings, in process technology, 
and in metal-production processes. In the geophysical and astrophysical context one 
may think of convection in Earth's mantle, in Earth's outer core, and in the outer layer 
of the Sun. E.g., the random reversals of Earth's or the Sun's magnetic field have been 
connected with thermal convection. 

Major progress in the understanding of the Rayleigh-Benard system has been 
made over the last decades, see e.g. the recent reviews [H [2]. Meanwhile it has 
been well established that the general heat transfer properties of the system, i. e. 
A/n = MuiJZa^Vr) and Tie = 7le{f/u,Vr), are well described by the Grossmann-Lohse 
(GL) theory [31 HI |5l |6] . In that theory, in order to estimate the thicknesses of the kinetic 
and thermal boundary layers (BL) and the viscous and thermal dissipation rates, the 
boundary layer flow is considered to be scalingwise laminar Prandtl-Blasius flow over a 
plate. We use the conventional definitions: The Rayleigh number is TZa = agH^A/uK 
with the isobaric thermal expansion coefficient a, the gravitational acceleration g, the 
height H of the RB system, the temperature difference A between the heated lower 
plate and the cooled upper plate, and the material constants u, kinematic viscosity, 
and K, thermal diffusivity, both considered to be constant in the container (Oberbeck- 
Boussinesq approximation). The Prandtl number is defined as Vr = u/k and the 
Reynolds number TZe = UH/u, with the wind amplitude U which forms in the bulk of 
the RB container. 

The assumption of a laminar boundary layer will break down if the shear 
Reynolds number TZcs in the BLs becomes larger than approximately 420 [7]. Most 
experiments and direct numerical simulations (DNS) currently available are in regimes 
where the boundary layers are expected to be still (scalingwise) laminar, see [1]. 
Indeed, experiments have confirmed that the boundary layers scalingwise behave as 
in laminar flow [8], i.e., follow the scaling predictions of the Prandtl-Blasius theory 
il [101 [III HSl HSl [3- Recently, Zhou et al. [H US] have shown that not only the scaling 
of the thickness, but also the experimental and numerical boundary layer profiles in 
Rayleigh-Benard convection agree perfectly with the Prandtl-Blasius profiles, if they 
are evaluated in the time dependent reference frames, based on the respective momentary 
thicknesses. This confirms that the Prandtl-Blasius boundary layer theory is indeed the 
relevant theory to describe the boundary layer dynamics in Rayleigh-Benard convection 
for not too large TZcs- 

The aim of this paper is to explore the consequences of the Prandtl-Blasius theory 
for the required numerical grid resolution of the BLs in DNSs. Hitherto, convergence 
checks can only be done a posteriori, by checking whether the Nusselt number does 
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not considerably change with increasing grid resolution [161 113 [IHl UHl EQl ET] or by 
guaranteeing (e.g. in ref. [22| [21]) that the Nusselt numbers calculated from the global 
energy dissipation rate or thermal dissipation rate well agree with that one calculated 
from the temperature gradient at the plates or the ones obtained from the overall 
heat flux. The knowledge that the profiles are of Prandtl-Blasius type offers the 
opportunity to a priori determine the number of required grid points in the BLs for 
given Rayleigh number and Prandtl number, valid in the boundary layer dominated 
ranges of moderately high TZa numbers. 

In section 2 we will first revisit the Prandtl-Blasius BL theory - see refs. 
[SI [ini im [121 1131 E] or for more recent discussions in the context of RB refs. [51 
- and derive the ratio between the thermal boundary layer thickness 6e and the velocity 
boundary layer thickness 6^ as functions of the Prandtl number Vr extending previous 
work (section 3). We will also discuss the limiting cases for large and small Vr, 
respectively. The transitional Prandtl number between the two limiting regimes turns 
out to be surprisingly small, namely Vr* = 0.046. The crossover range is found to 
be rather broad, roughly four orders of magnitude in Vr. In section 4 we note that 
the Prandtl-Blasius velocity BL thickness is different from the velocity BL thickness 
based on the position of the maximum r.m.s. velocity fluctuations (widely used in the 
literature), but well agrees with a BL thickness based on the position of the maximum 
of an energy dissipation derivate that was recently introduced in ref. [HI [2T]. We then 
derive the estimate for the minimum number of grid points that should be placed in 
the boundary layers close the top and bottom plates, in order to guarantee proper grid 
resolution. Remarkably, the number of grid points that must have a distance smaller 
than Su from the wall increases with increasing TZa, roughly as ~ TZa^'^^. This estimate 
is compared with a posteriori results for the required grid resolution obtained in various 
DNSs of the last three decades, finding good agreement. Section 5 is left to conclusions. 

2. Prandtl boundary layer equations 

The Prandtl-Blasius boundary layer equations for the velocity field u(x, z) (assumed to 
be two-dimensional and stationary) over a semi-infinite horizontal plate [9l [TOl [HI [121 



with the boundary conditions Ux{x,0) = 0, Uz{x,0) = 0, and Ux{x,oo) = U. Here 
Ux{x,z) is the horizontal component of the velocity (in the direction x of the large- 
scale circulation), Uz{x,z) the vertical component of the velocity (in the direction z 
perpendicular to the plate), and U the horizontal velocity outside the kinetic boundary 
layer (wind of turbulence). Correspondingly, the equation determining the (stationary) 
temperature field T(x, z) reads 



M, IZ] read 



UxdxUx + UzdzUx = udzdzUx, 



(1) 



UxdxT + Uzd.T = ndzdzT, 



(2) 
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with the boundary conditions T{x, 0) = Tpiate and T{x, oo) = Tbuik, which under 
Oberbeck-Boussinesq conditions is the arithmetic mean of the upper and lower plate 
temperature. Applying these equations to RB flow implies that we assume the 
temperature fleld to be passive. 

The dimensionless similarity variable ^ for the vertical distance z from the plate 
measured at the distance x from the plate's edge is 

e = zJ^. (3) 
V xu 

Since the flow in Prandtl theory is two-dimensional, a streamfunction \i/ can 
be introduced, which represents the velocity fleld. The streamfunction is non- 
dimensionalized as \1' = "^/y/xuU, and the temperature is measured in terms of A/2, 
giving the non-dimensional temperature fleld G. Rewriting eqs. ([T]) and (|2| in terms of 
and one obtains 

d^^/rf^^ + 0.5^rf2^/rf^2 =0, (4) 

d^e/d^'^ + 0.5Vr^de/d^ = 0. (5) 

Here the boundary conditions are 

^(0) = 0, d'^/d^iO) = 0, d<iJ/d^{oo) = 1, (6) 

0(0) = 0, e(oo) = l. (7) 

The temperature and velocity proflles obtained from numerically solving equations 
(|4])-([7]) (for particular Prandtl numbers) are already shown in textbooks [121 El US] and 
in the context of RB convection in refs. [231 [25]: From the momentum equation ^ with 
above boundary conditions one immediately obtains the horizontal velocity d'^ / d^. The 
dimensionless kinetic boundary layer thickness 6u can be deflned as that distance from 
the plate at which the tangent to the function d'^ / d^ at the plate (^ = 0) intersects the 
straight line d^! / d^ = 1 (see flgure[l]a). As equation Q and the boundary conditions ^ 
contain no parameter whatsoever, the dimensionless thickness 6u of the kinetic boundary 
layer with respect to the similarity variable ^ is universal, i.e., independent of Vr and 
U or Tie, 

5u = A-^ ^ 3.012 or A^ 0.332. (8) 

Solving numerically equation ^ with the boundary conditions ([T]) for any flxed 
Prandtl number, one obtains the temperature proflle with respect to the similarity 
variable ^ (see flgure[l]6). Note that in contrast to the longitudinal velocity d'^/d^, the 
temperature proflle G depends not only on ^ but also on the Prandtl number, since Vr 
appears in equation ([s]) as the (only) parameter. The distance from the plate at which 
the tangent to the proflle intersects the straight line B = 1 deflnes the dimensionless 
thickness of the thermal boundary layer, 

~5e = C{Vr), (9) 

where C{Vr) is a certain function of Prandtl number. E.g., one numerically flnds 
C ^ 3.417, 1.814, and 1.596 for Vr = 0.7, 4.38, and 6.4, respectively (see flgure |l| 
b). 
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Figure 1. Solution of the Prandtl-Blasius equations (|4])-([7|): (a) Longitudinal 
velocity profile ^(C) (solid curve) with respect to the similarity variable ^. The 
tangent to the longitudinal velocity profile at the plate = 0) and the straight line 
d^/d^ = 1 (both dashed lines) intersect at ^ = Su = A^^ w 3.012, for aU Vr. We define 
this value i5„ as the thickness of the kinetic boundary layer. (6) Temperature profile 
Q{£,) as function of the similarity variable ^ for Vr = 0.7 (black solid curve), Vr = 4.38 
(red solid curve) and Vr = 6.4 (green solid curve). The tangents to the profile curves 
at the plate (^ = 0) and the straight line 8 = 1 (dashed lines) define the edges 
(thicknesses) of the corresponding thermal boundary layers, i. e., ^ = 5^ = C{Vr). For 
the presented cases Vr = 0.7, 4.38, and 6.4 one has C(0.7) « 3.417, C(4.38) « 1.814, 
and C(6.4) « 1.596, respectively. 



From ([s]) and ^ one obtains the ratio between the (dimensional) thermal boundary 
layer thickness 5e and the (dimensional) kinetic boundary layer thickness 5u'- 

^i = Y = ^C{Vr). (10) 

As discussed above, the constant A and the function C = C{Vr) are found from the 
solutions of equations (|4])-([7]) for different Vr. A and C{Vr) reflect the slopes of the 
respective profiles, 



A 



(0), C{Vr) 



dQ 

~d^ 



1 -1 



(0) 



generally 



With (3) the physical thicknesses are 6u = ^u/ y ^ and 6e = Se/-\ 
depending on U and the position x along the plate. The physical thermal BL thickness 
then is 



CjVr) 

7W 



XV] 




U dQ 



-\ -1 



xu 



(0) 



dQ 

dz 



(0) 



(12) 



Thus, explicitly it depends neither on U nor on the position x along the plate. Reminding 
the definition of the thermal current J = {uzT) — Kdz{T), we get (^(0)) = ^72^^^^)) ~ 
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^ J = 2H ^Afu, i. e., on x-average we have 

Sg is the so-called slope thickness, see Sect. 2.4 of reference [23]. In contrast to the 

thermal BL thickness 6e the physical velocity BL thickness Su = I \J~^ depends 
explicitly both on the position x and on the wind amplitude U . In a Rayleigh-Benard 
cell we choose for x a representative value x = dL = dVH. Then the famous Prandtl 
formula [9] results 

^ (14) 



Tie 



Here a = ^^J ^ = A ^y/oT. The constant a has been obtained empirically based on 
the experimental measurements by [26] performed in a cylindrical cell of aspect ratio 
one, filled with water. The result was [5] 

a ^ 0.482. (15) 

We note that this value probably depends on the aspect ratio, on the shape of the 
RB container, and can also be different for numerical 2D Rayleigh-Benard convection 
[27] [25| I2S]. It will also be different for the slope thickness as considered here or other 
definitions as e. g. the 99% -thickness. 

It seems worthwhile to note that similarly to the case of dg also 5^ can be expressed 
by a profile slope at the plate. Analogously to the temperature case one calculates 
for the kinetic thickness 5^ = t//^(0). Here U appears explicitly and the derivative 
may depend on x. The denominator is the local stress tensor component, which - 
after averaging - describes the momentum transport, just as the temperature profile 



derivative at the plate characterises the heat transport. In combination with eq. (14) it 
says that the kinetic stress behaves as (%^(0)) ~ UVTZe/{aH). 



From eqs. (10) and (14) we also find the useful (and known) relation for Prandtl- 



Blasius boundary layers 



Sg = agC{Vr)^= with ag = A ■ a ^ 0.160. (16) 



From solving equations (|4])-([7]) together with relations (11) one obtains that the BL 
thickness ratio (10) has two limiting cases, namely Sg/Su ~ Vr~^^'^ for very small Pr ^ 1 
and 6g/6u ~ Vr~^^^ for very large Vr ^ 1. We thus present the ratio of the thermal 
and kinetic boundary layer thicknesses normalised by Vr^^^^ in figure [2] for different Vr 
from Vr = 10^^ to 10^. The figure confirms that the scaling of the ratio between the 
thermal and kinetic boundary layer thicknesses in the low and high Prandtl number 
regimes is Vr~^^'^ and Vr~^^^, respectively. Between these two limiting regimes there 
is a transition region, whose width is about 4 orders of magnitude in Vr. In the next 
section we will derive analytic expressions for the ratio 6g / 6u in the respective regimes, 
which will be used in the remainder of the paper to analyse the resolution properties 
for DNS in the BLs of the Rayleigh-Benard system. 
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Figure 2. Double- logarithmic plot of the ratio of the thermal and kinetic boundary 

1/3 

layer thicknesses, normalised by Vr ' , as obtained from numerical solution of 
equations Q-Q as function of Vr (solid black line). For large Vr the curve through 
the data is constant, for small Vr the (plotted, reduced) curve behaves cx Vr^^^^. 
Approximation ( 22 ) (green dotted line) is indistinguishable from Se /6u in the region 
Vr < 3x lO""*. Approximation (24) (blue dashed-dotted line) well represents Sg/5u for 
Vr > 0.3; for Vr > 3 it practically coincides with approximation (25 1. Approximation 



( 26 1 (red solid curve) connects the analytical approximations in the transition range 
3 X 10^^ < Vr < 3 between the lower and upper Prandtl number regimes. 



In the Prandtl-Blasius theory the asymptotic velocity amphtude f/ is a given 
parameter; the resulting heat current A/u is a performance of the boundary layers only. 
In contrast, in Rayleigh-Benard convection the heat transport is determined by the 
BLs together with the bulk flow. Therefore in RB convection the wind amplitude U no 
longer is a passive parameter, but U and J\fu are actively coupled properties of the full 
thermal convection process. 

The Reynolds number TZe is deflned as the dimensionless wind amplitude, 

7^e=^. (17) 



From the law for the kinetic BL thickness (14), the thermal BL thickness Sg (13), and 
the BL thickness ratio (10) one obtains 

This TZe ~ A/u^ law is in perfect agreement with the GL theory [3], HI El E]. In that 
theory several sub-regimes in the (Jla, Vr) parameter space are introduced, depending 
on the dominance of the BL or bulk contributions. In regimes I and II the BL of the 
temperature fleld dominates, while in III and VI it is the thermal bulk. Regimes I 
and II differ in the velocity fleld contributions: It either is the m-BL (I) or the w-bulk 
(II) which dominates; analogously the pair III and IV is characterized. The labels £ 
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(for lower Pr) and u (for upper Pr) distinguish the cases in which the thermal BL is 
thicker or smaller than the kinetic one. All ranges in the GL theory, which are thermal 
boundary layer dominated, show the TZe ~ A/u^ behaviour, namely Ii, lu, Hi, IIu- In 
the thermal bulk dominated ranges of RB convection the relation between TZe and Mu 
is different. In ///„ we have TZe ~ Mi^^^, in IVi it is TZe ~ A/u, and in IVu also 
TZe ~ Nu^^^ holds; but here the Prandtl-Blasius result ( 18 ) is not applicable, since the 
heat transport mainly depends on the heat transport properties of the bulk. In the 
range Joo, although boundary layer dominated, also a different relation {TZe ~ An^) 
holds; here the upper and the lower kinetic BLs fill the whole volume and therefore 
there is no free flow outside the BLs, in contrast to the Prandtl-Blasius assumption of 
an asymptotic velocity with the LSC amplitude U . 



3. Approximations for the ratio 80 /5u of the temperature and velocity 
boundary layer thicknesses 

In this section we will derive analytical approximations for the ratio 5e / 5u for the three 
regimes identified in the previous section, cf. figure [2j We start by discussing the low 
(Vr < 3 X 10~^) and the high (3 < Pr) Prandtl number regimes, before we discuss the 
transition region 3 x 10~^ < < 3. 



3.1. Approximation of 5e/5u forT'r < 3 x 10 



In the case of very small Prandtl number, Pr <^ 1, the thickness of the velocity boundary 
layer is negligible compared with the thickness of the temperature boundary layer, i.e., 
S> 5u- Hence, in most of the thermal boundary layer it is ~ U . Introducing the 
similarity variable as in ref. [13] 

_z [IT 

one obtains the following equation for the temperature as a function of rj: 



(19) 



d^Q/drj^ + 2ride/dri = 0, with 6(0) = 0, 6(00) = 1. 
The solution of this boundary value problem is the Gaussian error function 



6(77) = erf(?7) 



'dt. 



(20) 



According to (|3| and (19), the similarity variable ^ used in the Prandtl equations and 
the similarity variable 77 used in the approximation for Pr ^ 1 are related as follows 

= ipr^/2^. (21) 



Applying now the formulae (20), (21) and (11) we obtain the following equalities: 

2 de 



drj 



(0) 



— ((]) ^ 
dC dr] 



CiVr) 
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This leads to the approximation for the function C{Vr) = ^/^TVr~^^'^ for very small Vr. 

A^Vr-^l^ ^ 0.588Pr-^/^ Vr < 1. (22) 



In figure |2] one can see that for very small Prandtl numbers, Vr < 3 x 10~^, the 
approximation ([22|) is as expected indistinguishable from the numerically obtained 5e/Su- 



3. 2. Approximation of 60 / 6u for Vr > 3 

Meksyn [12j, based on the work by Pohlhausen [11], derived that the solution of the 
temperature equation (|5|, together with relation ^ equals 

® (71) ^ ^ lo^^' e-^W^^rft, F{t) = £ nq)dq. (23) 

The constant D can be found as usual from the boundary condition at infinity and was 
approximated in pTl [T2] for "Pr > 1 as follows 

^ 0.478Pr^/3 X . 1 1 161 

D = . . , c{Vr) ^ 1 + —— - — — ^ + 



c{Vr) ' ' A5Vr 405Vr^ 601425Pr^ 



From this and (23) one derives 



c{Vr) d{i/^y ' dC CiVr)' 

This connects c{Vr) and C{Vr) as follows 

\/2 

CiVr) ^ — — c(Vr)Vr-^/^ ^ 2.959 c(Pr)Pr"^/^ 
^ ' 0.478 ^ ^ V / ' 

resulting in the approximation 

-± = ACiVr) = EVr-^'^ciVr), E ^ A^— ^ 0.982. (24) 
6u 0.478 ^ ^ 

For Pr ^ 1, the function c{Vr) approaches 1, hence C{Vr) ^ 2.959Pr~^''^, implying 

J- = EVr-^^^, Vr > 1. (25) 
0,,, 



In figure 2 the approximation (24) is presented as a blue dash-dotted curve. For Vr > 3 



the function {5e/5u)Vr^^^ almost coincides with the constant E. 

3.3. Approximation of Sq/Su in the crossover range 3 x 10~^ <Vr<3 

As one can see in figure |2| the approximation (22) well represents 5e/5u in the region 



Vr < 3x 10^^, while (25) is a good approximation oiSe/6u for Vr > 3. An approximation 
of the ratio of the thermal and kinetic boundary layer thicknesses in the transition region 
3 X 10~^ < Pr < 3 is obtained by applying a least square fit to the numerical solutions 
of the Prandtl-Blasius equations (|4])-([7|). One finds: 



As seen in figure [2| this relation is a good fit of the full solution in the transition regime. 



_ P^-o.357+o.o22iogPr^ 3 X 10"^ < Pr < 3. (26) 
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3.4- Summary 

For the ratio 60/6^ of the thicknesses of the thermal and kinetic boundary layers, which 



depends strongly (and only) on Vr, we find according to (22), (25), and (26) 

Ay^Vr^^/^, A ^ 0.332, Vr < 3 x 10-^ 
P^-o.357+o.o22iogPr^ 3 X lO^^ < Pr < 3, (27) 

EVr-^'^, E ^ 0.982, Vr > 3. 



The crossover Prandtl number Vr* between the asymptotic behaviours, cf. first 



and last line of (27), is defined as the intersection point Vr* = 0.046 of the 
asymptotic approximations. Note that this crossover between the small-Pr behaviour 
6e/Su oc Vr^^^^ and the large- Pr behaviour 6e/6u oc Vr^^^^ does not happen at a Prandtl 
number of order 1, but at the more than 20 times smaller value Vr* = 0.046. In this 
sense most experiments are conducted in the large Vr regime. However, also note that 
other definitions of the BL thicknesses lead to other crossover Prandtl numbers. 

Finally, we also give the thickness of the kinetic BL in the three regimes, as obtained 



from (27) and (13), namely 

0.5Afu-^Vr^^^A-^n-^/^H, Vr < 3 x 10-^ 

0.5ArM-W-3"-°-°22iogPr^^ 3xlO-4<Pr<3, (2J 

Q.bNu'^Vr^/'^E'^H, Vr > 3. 



We compare this Prandtl-Blasius result (28) for the kinetic boundary layer 
thickness in terms of TVu and Vr (thus valid if the heat transport is BL dominated) 
with the estimate given in reference [21], where the kinetic boundary layer thickness in 
a cylindrical cell is identified as two times that height at which the averaged quantity 

e::=(u-V2u),,<^,, (29) 

has a maximum, because it was empirically found that the maximum of is 
approximately in the middle of the velocity boundary layer. Here u is the velocity 
field and the averaging is over time t, the azimuthal direction 0, and over the radial 
direction O.li? < r < 0.9-R, with R the radius of the cylindrical convective cell. The 
restricted range for the radial direction has been used in order to exclude the singularity 
region close to the cylinder axis and the region close to the sidewall, where the definition 
misrepresents the kinetic boundary layer thickness. Figure [3] shows that there is a very 
good agreement between the theoretical Prandtl-Blasius slope boundary layer thickness 



and that obtained using (29). The figure also shows that the position of the maximum 
r.m.s. velocity fluctuations is not a good indicator for the velocity boundary layer edge; 
it rather seems to identify the position where the LSC is the strongest. 



4. Resolution requirements within the boundary layers in DNS 

We now come to the main point of the paper: What can we learn from the Prandtl- 
Blasius theory for the required mesh resolution in the BLs of DNS of turbulent RB 
convection? Obviously, a "proper" mesh resolution should be used in order to obtain 
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Figure 3. Profiles of e„ (29) (black), and the r.m.s. velocity fluctuations for the 
azimuthal velocity component (green) for (a) TZa — 10® and Vr = 6.4 and (6) 
TZa = 2 X 10^ and Vr = 0.7. The profiles have been normalised with the respective 
maxima for clarity. The vertical black lines indicate the velocity boundary layer 



thickness based on (28 1. The red dashed and solid lines indicate the heights at which 
the quantity (29 1 has a maximum and two times this height, respectively. The 
vertical green line indicates the position of the maximum r.m.s. velocity fluctuations. 



accurate results. In a perfect DNS the local mesh size should be smaller than the local 
Kolmogorov r7/f(x, t) and Batchelor ^^^(x, t) scales (see e.g. ref. [30]), and the resolution 
in the boundary layers should be also sufficient, see e.g. [311 [161 ISSl |32l [21]. It indeed 
has been well established that the Nusselt number is very sensitive to the grid resolution 
used in the boundary layers; when DNS is underresolved, the measured Nusselt number 
is too high [nU [331 HSl [311 133 [3S1 EI]- Hitherto, the standard way to empirically check 
whether the mesh resolution is sufficient is to try a finer mesh and to make sure that the 
Nusselt number is not too different. In this way the minimal number of grid points that 
is needed in the boundary layer is obtained by trial and error: Grotzbach [31j varied 
the number of grid points in the boundary layer between 1 and 5 in simulations up to 
TZa = 3 X 10^ with Vr = 0.71 and found that 3 grid points in the boundary layers should 
be sufficient. Verzicco and Camussi [33] tested this at TZa = 2 x 10^ and Vr = 0.7 and 
stated that at least 5 points should be placed in the boundary layers. Stevens et al. 
[21] tested the grid resolution for T^a = 2 x 10^ to 2 x 10^^ and Vr = 0.7. They found 
that for 7?a = 2 X 10^ the minimum number of nodes in the boundary layers should be 
around 10 and that this number increases for increasing TZa. Together with the earlier 
series of papers the data clearly suggest that indeed there is an increase of required grid 
points in the BL with increasing Rayleigh number. 

However, one must be careful. The empirical determination of the required number 
of grid points in the BL is not only intensive in computational cost, but also difficult. 
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The Nusselt number obtained in the simulations not only depends on the grid resolution 
in the BLs at the top and bottom plates, but also on the grid resolution in the bulk and 
at the side walls where the thermal plumes pass along [21]. So obviously a general 
theory-based criterion for the required grid resolution in the thermal and kinematic 
boundary layers will be helpful for performing future simulations. In this section we 
will derive such a universal criterion, harvesting above results from the Prandtl-Blasius 
boundary layer theory. 

We first define the (local) kinetic energy dissipation rates per mass, 

z/ >^ >^ / dui{x,t) duj{x,t) 

Its time and space average for incompressible flow with zero velocity b.c. is {eu)t,v = 

^ Xli Ylj ( ^'^dx'^^ ) )ty- connected with the Nusselt number through the exact 

relation 




{eu)t.y = —^{Mu- l)'RaVr-\ (31) 

This follows directly from the momentum equation for Rayleigh-Benard convection in 
Boussinesq approximation |37]. Here, {■)t,v denotes averaging over the whole volume of 
the convective cell and over time and (later) {■)t,A denotes averaging over any horizontal 
plane and time. 

We start with the well established criterion that in a perfect DNS simulation the 
(local) mesh size must not be larger than the (local) Kolmogorov scale [M] '7_ft'(x, t), 
which is locally defined with the energy dissipation rate of the velocity, 

r/^(x,t)= (z.V6„(x,t))'/\ (32) 

TjK is the length scale at which the inertial term u1/r and the viscous term ~ vu^jr'^ 
of the Navier-Stokes equation balance, where Ur ~ (e^?")^'''^ has been assumed for the 
velocity difference at scale r. A corresponding length scale t^t follows from the balance 
of the advection term ~ UrT^jr and the thermal diffusion term uTr/r"^ in the advection 
equation; it is 

r/T(x,t) = (/cVe„(x,t))'/' = r/^(x,t)Pr-3/^ (33) 

However, for large Pr the velocity field is smooth at those scales at which the temperature 
field is still fluctuating. Then the velocity difference Ur ~ \/^u/^f and advection term 
and thermal diffusion term balance at the so-called Batchelor scale [39] t^b, which is 
defined as 

r/B(x,t) = (z//«Ve„(x,t))'/' = r/^(x,t)Pr-^/2. (34) 

For small Vr < 1 obviously tjt > tjb > tjk and for comparison with the grid resolution, 
the Kolmogorov scale tjk seems to be the most restrictive (i.e., smallest) length scale. 
In contrast, for large Vr > 1 it rij- < tjb < rjx and one may argue that tjt is the most 
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restrictive length scale. This indeed may be the case in the Prandtl number regime in 
which the velocity field can still be described through Kolmogorov scaling Ur ~ (e„r)^/^, 
but for even larger Vr the velocity field becomes smooth Ur ~ \/ ^ul^r and then the 
grid resolution should be compared to the Batchelor scale r]B as smallest relevant length 
scale. In below analysis, for Pr > 1 we will restrict ourselves to this limiting case. 



We now define global Kolmogorov and Batchelor length scales 7^^°''"' = 



global ^ 1/1/4^1/2 
'IB - , a/4 ) 



respectively, (and also the global length scale 7]^°'^"'^ 



y3/4 

/, \l/4 



and 



Using 



the exact relation (31), one can find how the global Kolmogorov length 7^^°*"' depends 
on TZa, Vr, and Afu, namely 

^3/4 p^l/2 

(35) 



global 

Vk 



(e )^'^ 



-H. 



The admissible global mesh size h^^"^"''- should clearly be smaller than both -q^^^""^ and 
^gUibai ^ which implies that one is on the safe side provided that 

p^l/2 



obal 



< 



na^'\j\fu - lyi' 



H for Vr <l 



or with the relation (34) between the Kolmogorov and Batchelor length 

1 



obal 



< 



-H for Vr > 1. 



(36) 



(37) 



A similar way to estimate mesh requirements in the bulk was suggested for the first 
time by Grotzbach |31] . Note that with these estimates for the required bulk resolution 



for most times and locations one is on the safe side, as equation (31) is an estimate for 
the volume averaged energy dissipation rate, which is localized in the boundary layers. 
However, not only the background field but also plumes detaching from the boundary 
layers do require an adequate resolution. 

To estimate the number of nodes that should be placed in the boundary layers, we 
will first estimate the area averaged energy dissipation rate in a horizontal plane in the 
velocity BL, {eu)t,AeBL- Employing eqs. (17), (14) and (30), one can find a lower bound 
for this quantity, namely 

2^ 



A&BL 



> V 




> V 



dz 



t,A 



V 



t,A, 



(38) 



From eqs. (31), (38), (18) and (27) it follows a lower bound for the ratio 

6 



{eu)t,A&BL_ ^ Vr 7^e 



(en) 



ty 



o^VaNu 



64a 



,Vr' 
TZa 



647c^a^A'^Afu^Vr-^TZa-\ 
64aWM^Pr-°-i5+o-i=^2iogPr^- 

64a^E^Aru^TZa~\ 



Vr <3x 10-^ 
3 X 10-^ <Vr <3, 
Vr > 3. 



(39) 
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For the Kolmogorov length rj^ in the velocity BL one can therefore write 




1/4 



global 

Vk ■ 



(40) 

t,A<^BL 

The mesh size hP^ in the BL must be smaller than and ?7^^, i.e., one is on the safe 
side if 



2-3/2a-iA/5i-3/2pr3/4A-3/27r-3/4i7, < 3 X 10-^ 

2-3/2a-iA/5i-3/2prO-5355-o.033iogPr^^ 3 X 10-4 < Pr < 1, 
2-3/2a-W?i-'/2prO-0355-o.033iogPr^^ 1< < 3, 



(41) 



according to (l39|, (|40j), (jSej) and (|37j). 

From the relations (41), (27) and (13) one can estimate the minimum number of 
nodes of the computational mesh, which must be placed in each thermal and kinetic 
boundary layer close the plates. We find that this minimum number of nodes in the 
thermal boundary layers is 

< 3 X 10"^ 



> 



(42) 



> 



(43) 



v^aArni/2p^-o.5355+o.033iogPr^ 3 ^ ^^-4 <pr<l, 
v^aA/-ni/2p^-o.o355+o.o33iogPr^ 1< Pr < 3, 
^ V^aUv'/^E^I^ Vr>3, 
while the minimum number of nodes in the kinetic boundary layers is 

= - 

V.BL - f^BL - S^hBL 

v^aArMV2p^-o.i785+o.oiiiogPr^ 3 ^ ^q-a <Vr <l, 
y2aArui/2p^o.32i5+o.oiiiogPr^ 1< Pr < 3, 
^ v^aAT^i/^pr 1/3^1/2, Pr>3. 

The number of nodes in the thermal boundary layer looks very restrictive for very low 
"Pr; however, one should realise that for very low Vr the thermal boundary layer also 
becomes much thicker than the velocity boundary layer. Hence, the criterion for the 
number of nodes in the thermal boundary layers determines the ideal distribution of 
nodes above the viscous boundary layer. For very high Vr the kinetic boundary layer 
becomes much thicker than the thermal boundary layer, and hence the restriction for 
the velocity boundary layer determines the ideal distribution of nodes above the thermal 
boundary boundary layer. Note that for large Vr equation (42 ) suggests that the number 
of grid points in the thermal boundary layer becomes independent of Vr (for fixed Mu): 
Indeed, as the velocity field is smooth anyhow, with increasing Vr no extra grid points 
are necessary in the thermal BL. 

In figure|4]we show the minimum number of nodes A^th . BL ^^"^ . BL' respectively, 
necessary to simulate the cases which have been investigated experimentally so far, for 



Boundary layer structure in turbulent thermal convection 



15 



(a) 
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Figure 4. Minimum number of BL nodes necessary in DNS of boundary layer 



dominated, moderately high RB convection, (a) -^th.BL (42) in the thermal boundary 



layers and (&) _-qL (43) in the kinetic boundary layers, required to simulate the 
experimentally investigated cases, references 00] (lilac squares, Vr = 0.67), [21] (black 
triangles, 0.60 <Vr < 7.00), [22] (blue circles, 0.68 <Vr < 5.92), [23] (green triangles, 
0.73 < Vr < 6.00), [22 (red pentagons, 3.76 < Vr < 5.54), [45] (black crosses, 
Vr — 4.2) and [46] (black pluses, Vr — 7.0). Dashed lines are fits to the quasi-data 
(measured values introduced into eqs. (42 1, (43)), with precision 0(10^''); rounding 
the respective numbers to their upper bounds gives 

(a) ^th.BL « 0.357^a"■^^ (jiij) and 

0.67 <Vr< 0.73 



0.357^a" 

(liHl for the quasi-data in the ranges 10*^ < 72a < 10^" and 



different TZa and Vr. The data points are generated by introducing the experimental 
values of 7ia,Vr, and the (measured) corresponding Afu into the formulas (42), (43). 
Based on these quasi-data points, one can give e. g. the following fits for the minimum 
number of nodes within the boundary layers for the case of Vr ^0.7: 

^th.BL ~ 0.357?a°•^^ 10^ < Tea < 10^°, (44) 
^v.BL ~0•31^°•^^ 10^ < T^a < 10^°. (45) 
Note that the numerical pre-factors in these estimates significantly depend on the 
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Prandtl number and on the empirically determined (ref. [5]) value of a, cf. eq. (15) 



The minimum node numbers for other values of Vr can be calculated directly using the 



relations (42)-(43). Apparently the scaling exponent depends much less on Vr. - All 
these estimates only give lower bounds on the required number of nodes in the boundary 
layers. 

As discussed at the beginning of this section, previous studies by Grotzbach [31] . 
Verzicco and Camussi [33] , and Stevens et al. [21j found an increasing number of nodes 
that should be placed in the thermal and kinetic boundary layers. The theoretical results 
thus confirm all above studies, because the increasing number of nodes was due to the 
increasing TZa number at which the tests were performed. To be more specific: according 



to the estimates (44) and (45) for Vr = 0.7 the minimum number of nodes that should 
be placed in the thermal and kinetic boundary layers is ~ 2.3 for TZa = 3 x 10^, 
A^ fa 4.4 for TZa = 2 x 10^, and A^ ^ 8.7 for TZa = 2 x 10^. The empirically found values 
at the respective TZa with Vr ^ 0.7 are 3 for TZa = 3 x 10^, 5 for TZa = 2 x 10'', and 
10 for TZa = 2 x 10^. Thus there is very good agreement between the theoretical results 
and the empirically obtained values, especially if one considers the difficulties involved 



in determining these values empirically, and the empirical value for the constant a ( 15 ) 
that is used in the theoretical estimates. We want to emphasize that not only the 
boundary layers close to the plates, but also the kinetic boundary layers close to the 
vertical walls must be well resolved. 

To sum up, the mesh resolution should be analysed a priori using the resolution 



requirements in the bulk (36), (37) and in the boundary layers (42), (43). Having 
conducted the DNS, the Kolmogorov and Batchelor scale should be checked a posteriori, 
to make sure that the mesh size was indeed small enough (as it has been done, for 
example, in refs. [T9l 120] ). 

5. Conclusion 

In summary, we used laminar Prandtl-Blasius boundary layer theory to determine the 



relative thicknesses of the thermal and kinetic boundary layers as functions of Pr (27). 

We found that neither the position of the maximum r.m.s. velocity fluctuations 
nor the position of the horizontal velocity maximum reflect the slope velocity boundary 
layer thickness, although many studies use these as criteria to determine the boundary 
layer thickness. In contrast to them, the algorithm by Stevens et al. [21] agrees very 
well with the theoretical estimate of the kinetic slope boundary layer thickness. 

We used the results obtained from the Prandtl-Blasius boundary layer theory to 
derive a lower bound on the minimum number of nodes that should be placed in the 
thermal and kinetic boundary layers close to the plates. We found that this minimum 
number of nodes increases not slower than ~ TZa^'^^ with increasing TZa. This result is in 
excellent agreement with results from several numerical studies over the last decades, in 
which this minimum number of nodes was determined empirically. Hence, the derived 
estimates can be used as guideline for future direct numerical simulations. 
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